#include "BWT.cu"
#include "cuda_header.h"
#include "mian_util.h"
#include "compaction.cu"
#include <thrust/scan.h>
#include <thrust/reduce.h>
#include "thrust/sort.h"
#include "thrust/device_vector.h"
#include "thrust/host_vector.h"
#include <cuda.h>
//#include"OccValue_Lam.cu"

__device__ int BackwardSearch1Error2(uint3* result, int maxNumOutput,
		BWT * bwt, unsigned int* theQuery,
		unsigned int queryId,
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_start, unsigned int * sa_end, 
		long* numAlign = NULL)
{
	unsigned int temp_start=(*sa_start),temp_end=*sa_end;
	int i;
	char c, misBase;
	int numResult = 0;

	for (i=upper-1; i >=(int)lower && *sa_start<=*sa_end; i--)
	{
		c = getCharFromQueries(theQuery, i);

		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{		
			if (misBase == c)
				continue;
			temp_start = *sa_start;
			temp_end = *sa_end;

			if(!bwtCalc(bwt,misBase,&temp_start,&temp_end))
				continue;

			//add pausedQuery to buffer	
			if(BackwardSearch(bwt, theQuery, 0,i, &temp_start, &temp_end))
			{
				if(numResult + 1 < maxNumOutput) {

					//if(temp_end - temp_start + 1 < 100) {
						result[numResult].x = queryId;
						result[numResult].y = temp_start;
						result[numResult].z = temp_end;
						numResult++;
					//}
				
					if(numAlign != NULL) {
						*numAlign += (temp_end - temp_start + 1);
						if(*numAlign > MAX_NUM_ALIGN_PER_QUERY) {
							*numAlign = -1;
							return 0;
						}
					}
				} else {
					return -1;
				}
			}
		}
		bwtCalc(bwt, c,sa_start,sa_end);
	}

	return numResult;
}



__device__ int ForwardSearch1Error2(uint3* result, int maxNumOutput,
		BWT * revbwt, unsigned int * query, unsigned int queryId,
		unsigned int lower, unsigned int upper, unsigned int mismatchCount, 
		unsigned int * sa_rev_start, unsigned int * sa_rev_end, 
		unsigned int * sa_start, unsigned int * sa_end, 
		unsigned int * sharedOccCount, long * numAlign = NULL)
{
	unsigned int temp_start=*sa_start,temp_end=*sa_end, temp_rev_start=*sa_rev_start, temp_rev_end=*sa_rev_end;
#if (OCC_COUNT_PER_THREAD >= 8)
	unsigned int *occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	unsigned int i;
	char c, misBase;
	int numResult = 0;


	for (i=lower; i<upper && *sa_start<= *sa_end; i++) {

		bwtCalcAll(revbwt,sa_rev_start, sa_rev_end,occCount_start, occCount_end, occCount);
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;

			temp_rev_start = *sa_rev_start;
			temp_rev_end = *sa_rev_end;

			bwtCalcWithOccValue(revbwt->countValue, misBase,occCount_start[misBase], occCount_end[misBase], 
					&temp_rev_start, &temp_rev_end);

			temp_end = *sa_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start > temp_end)
				continue;


			if(ForwardSearch(revbwt, query, i+1,upper, &temp_rev_start, &temp_rev_end, &temp_start,&temp_end,
						sharedOccCount + ALPHABET_SIZE))
			{
				if(numResult + 1 < maxNumOutput) {

					//if(temp_end - temp_start + 1 < 100) {
						result[numResult].x = queryId;
						result[numResult].y = temp_start;
						result[numResult].z = temp_end;
						numResult++;
					//}

					if(numAlign != NULL) {
						*numAlign += (temp_end - temp_start + 1);
						if(*numAlign > MAX_NUM_ALIGN_PER_QUERY) {
							*numAlign = -1;
							return 0;
						}
					}

				} else {
					return -1;
				}

				//ReportResult(mismatchCount+1, temp_start, temp_end, saInterval);
			}

		}

		bwtCalcWithOccValue(revbwt->countValue, c,occCount_start[c], occCount_end[c], sa_rev_start, sa_rev_end);

		*sa_end = *sa_end - occCount[c];
		*sa_start = *sa_end - (*sa_rev_end - *sa_rev_start);
	}

	return numResult;
}
__device__ int BWT1MismatchSearch_Backward(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId,
		BWT *bwt, BWT *revbwt, long * numAlign = NULL)
{
	unsigned int start= 1, end = bwt->textLength-1;
	char c;
	unsigned int i;
	int tempN = 0, numOutput=0;
	unsigned long long packedPattern = 0;
	int queryLen = QUERY_LENGTH;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, queryLen - TABLE_SIZE +i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];

	//      BackwardSearch(bwt, Queries, queryId, sizeX, QUERY_LENGTH-1, &start, &end);
	if(start <= end && BackwardSearch(bwt, query, (QUERY_LENGTH / 2), queryLen-TABLE_SIZE, &start, &end))
	{
		tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
				bwt, query, queryId,0,(QUERY_LENGTH / 2),0,
				&start, &end, numAlign);

		if( numAlign!=NULL &&  *numAlign == -1)
			return 0;
		if(tempN == -1) {
			return -1;
		} else {
			numOutput += tempN;
		}
	}
	return numOutput;
}


__device__ int BWT1MismatchSearch_Forward(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId,
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount, long * numAlign = NULL)
{

	unsigned int start= 1, end = bwt->textLength-1, revStart, revEnd;
	char c;
	unsigned int i;

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	int numOutput = 0;
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= ( c & 3 );
	}
	r_packedPattern = l_packedPattern;
	l_packedPattern <<= (TABLE_SIZE - revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (TABLE_SIZE -revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(TABLE_SIZE - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	if (end - start < revEnd - revStart) {
		revStart =revEnd - (end - start );
	} else if (end - start > revEnd - revStart) {
		start = end - (revEnd - revStart);
	}

	ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize, (QUERY_LENGTH / 2), &revStart, &revEnd, &start, &end,
			sharedOccCount);

	if(start <= end){
		int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
				revbwt, query, queryId, (QUERY_LENGTH/2),QUERY_LENGTH,0,
				&revStart,&revEnd,&start, &end, sharedOccCount, numAlign);
		if(numAlign != NULL && (*numAlign) == -1)
			return 0;
		if(tempN == -1) {
			return -1;
		} else {
			numOutput += tempN;
		}
	}
	return numOutput;
}

__global__ void BWT1MismatchBackward_Kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;

		numOutputCurCase = BWT1MismatchSearch_Backward(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}

__global__ void BWT1MismatchBackward_limit_Kernel(long* d_numAlign, uint3* d_result, int* d_numOutput,
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurCase;
	BWT * ptr_bwt = d_bwt;
	BWT * ptr_revbwt = d_revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		if(d_numAlign[idx] == -1 || d_numOutput[idx] == -1)
			continue;	
		numOutputCurCase = BWT1MismatchSearch_Backward(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt, d_numAlign + idx);

		if(d_numAlign[idx] == -1) {
			continue;
		}
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} 
		
		d_numOutput[idx] += numOutputCurCase;
	}
}

__global__ void BWT1MismatchForward_Kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{

	__shared__ unsigned int sharedOccCount[MAX_THREAD_NUM * OCC_COUNT_PER_THREAD];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;

		numOutputCurCase = BWT1MismatchSearch_Forward(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt, 
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}
__global__ void BWT1MismatchForward_limit_Kernel(long* d_numAlign, uint3* d_result, int* d_numOutput,
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	__shared__ unsigned int sharedOccCount[MAX_THREAD_NUM * OCC_COUNT_PER_THREAD];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurCase;
	BWT * ptr_bwt = d_bwt;
	BWT * ptr_revbwt = d_revbwt;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		if(d_numAlign[idx] == -1 || d_numOutput[idx] == -1)
			continue;
		numOutputCurCase = BWT1MismatchSearch_Forward(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD, d_numAlign + idx);

		if(d_numAlign[idx] == -1)
			continue;
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} 
		
		d_numOutput[idx] += numOutputCurCase;
	}
}

__device__ int BWT2MismatchSearch_CaseA2(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt)
{
	//unsigned int* query = Queries + QUERY_SIZE_IN_WORD*queryId;
	unsigned int temp_start=0,temp_end=0;

	//unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = floor((double)QUERY_LENGTH * 0.3f) * 2;

	unsigned long long packedPattern = 0;

	unsigned int i;
	char c, misBase;
	unsigned int start, end;
	//	unsigned int mtableSize = bwt->lookupTable->tableSize;
	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, QUERY_LENGTH - TABLE_SIZE + i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];


	int numOutput = 0;
	int tempN = 0;

	if(start <= end &&  BackwardSearch(bwt, query, sizeY, QUERY_LENGTH - TABLE_SIZE, &start, &end))
	{
		for (i=sizeY-1; i>=1 && start<=end; i--)
		{
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if( misBase == c)
					continue;
				temp_start = start;
				temp_end = end;
				if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){
					tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
							bwt, query, queryId,0,i,1, &temp_start, &temp_end);
					if(tempN == -1) {
						return -1;
					} else {
						numOutput += tempN;
					}
				}
			}
			bwtCalc(bwt, c, &start, &end);
		}

	}

	return numOutput;
}



__device__ int BWT2MismatchSearch_CaseA2_limit(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, long* numAlign)
{
	unsigned int temp_start=0,temp_end=0;
	unsigned long long packedPattern = 0;

	unsigned int i;
	char c, misBase;
	unsigned int start, end;
	for (i = 0; i < TABLE_SIZE ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, QUERY_LENGTH - TABLE_SIZE + i);
		packedPattern |= (c & 3);
	}
	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];


	int numOutput = 0;
	int tempN = 0;

	if(start <= end &&  BackwardSearch(bwt, query, sizeY, QUERY_LENGTH - TABLE_SIZE, &start, &end))
	{
		for (i=sizeY-1; i>=1 && start<=end; i--)
		{
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if( misBase == c)
					continue;
				temp_start = start;
				temp_end = end;
				if(bwtCalc(bwt,misBase,&temp_start,&temp_end)){

					long numAlignNew = *numAlign;
					tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
							bwt, query, queryId,0,i,1, &temp_start, &temp_end, &numAlignNew);

					if(numAlignNew == -1) {//this query will be discarded
						*numAlign = -1;
						return 0;
					}

					if(tempN == -1) {
						return -1;
					} else {
						numOutput += tempN;
						*numAlign = numAlignNew;
					}
				}
			}
			bwtCalc(bwt, c, &start, &end);
		}
	}

	return numOutput;
}




__device__ int BWT2MismatchSearch_CaseB2(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount)
{
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int *occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long l_packedPattern = 0;
	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	//unsigned long long mask;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;
	int numOutput = 0;

	l_packedPattern = 0;
	r_packedPattern = 0;
	rev_packedPattern = 0;
	// look-up the first characters
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= (c & 3);
	}
	r_packedPattern = l_packedPattern;
	// if the look-up tables are of different size
	l_packedPattern <<= (bwt->lookupTable->tableSize -  revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end  = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end - start<revEnd-revStart)
	{
		revStart = revEnd - (end - start);
	} else if (end - start > revEnd - revStart)
	{
		start = end - (revEnd - revStart);
	}

	if (start <= end && ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize,sizeY,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE))
	{
		for (i=sizeY; i<sizeZ && start<=end; i++)
		{
			bwtCalcAll(revbwt, &revStart, &revEnd, occCount_start, occCount_end, occCount);
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if(misBase == c)
					continue;
				temp_rev_start = revStart;
				temp_rev_end = revEnd;

				bwtCalcWithOccValue(revbwt->countValue,misBase,
						occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

				temp_end = end - occCount[misBase];
				temp_start = temp_end - (temp_rev_end - temp_rev_start);

				if(temp_start <= temp_end){

					int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
							revbwt, query, queryId, i+1,sizeZ, 1, 
							&temp_rev_start, &temp_rev_end,&temp_start, &temp_end,
							sharedOccCount + ALPHABET_SIZE);
					if(tempN == -1) {
						return -1;
					}
					else {
						numOutput += tempN;
					}

				}
			}

			bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

			end = end - occCount[c];
			start = end - (revEnd - revStart);
		}
	}

	return numOutput;
}




__device__ int BWT2MismatchSearch_CaseB2_limit(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount, long *numAlign)
{
	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int*occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif
	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	unsigned long long l_packedPattern = 0;
	//unsigned long long mask;
	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;
	int numOutput = 0;

	l_packedPattern = 0;
	r_packedPattern = 0;
	rev_packedPattern = 0;
	// look-up the first characters
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		l_packedPattern<<=2;
		c = getCharFromQueries(query, i);
		l_packedPattern |= (c & 3);
	}
	r_packedPattern = l_packedPattern;
	// if the look-up tables are of different size
	l_packedPattern <<= (bwt->lookupTable->tableSize -  revbwt->lookupTable->tableSize)*2;
	r_packedPattern <<= (bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2;
	r_packedPattern |= (1<<(bwt->lookupTable->tableSize - revbwt->lookupTable->tableSize)*2) - 1;

	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	start = l_packedPattern ? bwt->lookupTable->table[l_packedPattern-1]+1 : 1;
	end  = bwt->lookupTable->table[r_packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end - start<revEnd-revStart)
	{
		revStart = revEnd - (end - start);
	} else if (end - start > revEnd - revStart)
	{
		start = end - (revEnd - revStart);
	}

	if (start <= end && ForwardSearch(revbwt, query, revbwt->lookupTable->tableSize,sizeY,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE))
	{
		for (i=sizeY; i<sizeZ && start<=end; i++)
		{	
			bwtCalcAll(revbwt, &revStart, &revEnd, occCount_start, occCount_end, occCount);
			c = getCharFromQueries(query, i);
			for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
			{
				if(misBase == c)
					continue;
				temp_rev_start = revStart;
				temp_rev_end = revEnd;

				bwtCalcWithOccValue(revbwt->countValue,misBase,
						occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

				temp_end = end - occCount[misBase];
				temp_start = temp_end - (temp_rev_end - temp_rev_start);

				if(temp_start <= temp_end){


					long numAlignNew = *numAlign;
					int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
							revbwt, query, queryId, i+1,sizeZ, 1, 
							&temp_rev_start, &temp_rev_end,&temp_start, &temp_end, 
							sharedOccCount + ALPHABET_SIZE, &numAlignNew);
					if(numAlignNew == -1){
						*numAlign = -1;
						return 0;
					}

					if(tempN == -1) {
						return -1;
					} else {
						numOutput += tempN;
						*numAlign = numAlignNew;
					}

				}
			}

			bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

			end = end - occCount[c];
			start = end - (revEnd - revStart);
		}
	}

	return numOutput;
}





__device__ int BWT2MismatchSearch_CaseC2(uint3* result, int maxNumOutput,
		unsigned int * query, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount)
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int*occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif

	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	rev_packedPattern = 0;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(query, i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(query, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3 );
	}

	//for error happens in lookup range
	for(i=sizeX; i<bwt->lookupTable->tableSize; i++)
	{
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if(misBase == c)
				continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;

			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start);
			}else if ( end - start > revEnd - revStart) {
				revStart = revEnd -( revEnd - revStart);
			}

			if(start > end)
				continue;
			unsigned int pos = revbwt->lookupTable->tableSize;
			ForwardSearch(revbwt, query, pos, sizeY, &revStart, &revEnd, 
					&start, &end, sharedOccCount + ALPHABET_SIZE);

			if( start <= end ){

				int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						revbwt,  query, queryId,sizeY, sizeZ, 1, 
						&revStart,&revEnd,&start, &end, sharedOccCount + ALPHABET_SIZE);

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
				}
			}
		}
	}




	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];





	//Lookup SA correction...
	if (end - start < revEnd - revStart) {
		revStart = revEnd - ( end - start);
	} else if (end - start > revEnd - revStart) {
		start = end - ( revEnd - revStart);
	}

	if(sizeX > revbwt->lookupTable->tableSize){
		ForwardSearch(revbwt, query, TABLE_SIZE ,sizeX,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE);
		i = sizeX;
	}
	else 
		i = TABLE_SIZE;
	for (; i<sizeY && start<=end ; i++)
	{
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		c = getCharFromQueries(query, i);
		for(misBase=0; misBase < ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if(temp_start > temp_end)
				continue;
			if(ForwardSearch(revbwt, query, i+1,sizeY,&temp_rev_start,&temp_rev_end,
						&temp_start,&temp_end, sharedOccCount + ALPHABET_SIZE) 
					&& temp_start <= temp_end) {

				int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						revbwt, query, queryId,sizeY,sizeZ,1,&temp_rev_start,
						&temp_rev_end,&temp_start, &temp_end, sharedOccCount + ALPHABET_SIZE);

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
				}
			}
		}

		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}

	return numOutput;
}




__device__ int BWT2MismatchSearch_CaseC2_limit(uint3* result, int maxNumOutput,
		unsigned int * Queries, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount, long* numAlign )
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int*occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif


	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long packedPattern = 0;
	unsigned long long rev_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	rev_packedPattern = 0;

	for (i = 0; i <bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(Queries, i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(Queries, revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3 );
	}

	//for error happens in lookup range
	for(i=sizeX; i<bwt->lookupTable->tableSize; i++)
	{
		c = getCharFromQueries(Queries, i);
		for(misBase=0; misBase<ALPHABET_SIZE; misBase++)
		{
			if(misBase == c)
				continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;

			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start);
			}else if ( end - start > revEnd - revStart) {
				revStart = revEnd -( revEnd - revStart);
			}

			if(start > end)
				continue;
			unsigned int pos = revbwt->lookupTable->tableSize;
			ForwardSearch(revbwt,Queries, pos, sizeY, &revStart, &revEnd, 
					&start, &end, sharedOccCount + ALPHABET_SIZE);

			if( start <= end ){

				long numAlignNew = *numAlign;

				int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						revbwt,  Queries, queryId,sizeY, sizeZ, 1, 
						&revStart,&revEnd,&start, &end, 
						sharedOccCount + ALPHABET_SIZE, &numAlignNew);
				if(numAlignNew == -1){
					*numAlign = -1;
					return 0;
				}

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
					*numAlign = numAlignNew;
				}

			}
		}
	}




	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];





	//Lookup SA correction...
	if (end - start < revEnd - revStart) {
		revStart = revEnd - ( end - start);
	} else if (end - start > revEnd - revStart) {
		start = end - ( revEnd - revStart);
	}

	if(sizeX > revbwt->lookupTable->tableSize){
		ForwardSearch(revbwt, Queries, TABLE_SIZE ,sizeX,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE);
		i = sizeX;
	}
	else 
		i = TABLE_SIZE;
	for (; i<sizeY && start<=end ; i++)
	{
		c = getCharFromQueries(Queries, i);
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		for(misBase=0; misBase < ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if(temp_start > temp_end)
				continue;
			if(ForwardSearch(revbwt, Queries, i+1,sizeY,&temp_rev_start,&temp_rev_end,
						&temp_start,&temp_end, sharedOccCount + ALPHABET_SIZE) 
					&& temp_start <= temp_end) {

				long numAlignNew = *numAlign;
				int tempN = ForwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						revbwt, Queries, queryId,sizeY,sizeZ,1,&temp_rev_start,
						&temp_rev_end,&temp_start, &temp_end,  
						sharedOccCount + ALPHABET_SIZE, &numAlignNew);


				if(numAlignNew == -1){
					*numAlign = -1;
					return 0;
				}

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
					*numAlign = numAlignNew;
				}
			}
		}

		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}

	return numOutput;
}





__device__ int BWT2MismatchSearch_CaseD2_limit(uint3* result, int maxNumOutput,
		unsigned int * Queries, unsigned int queryId, 
		BWT *bwt, BWT *revbwt, unsigned int * sharedOccCount,long* numAlign )
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int*occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif


	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];

	unsigned long long packedPattern = 0;
	//unsigned long long l_packedPattern = 0;
	//unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	//unsigned long long rev_l_packedPattern = 0;
	//unsigned long long rev_r_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	//l_packedPattern = 0;
	//r_packedPattern = 0;
	rev_packedPattern = 0;
	//rev_l_packedPattern = 0;
	//rev_r_packedPattern = 0;

	for (i = 0; i < bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	//For error happen in lookup range....
	for (i = sizeY - sizeX; i< bwt->lookupTable->tableSize ;i++) {
		c = getCharFromQueries(Queries, sizeX + i);
		for (misBase=0;misBase<4;misBase++) {
			if (misBase == c) continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;
			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start );
			} else if ( end - start >revEnd - revStart) {
				start = end - (revEnd - revStart);
			}
			if( start <= end && ForwardSearch(revbwt, Queries, sizeX+revbwt->lookupTable->tableSize, QUERY_LENGTH, 
						&revStart, &revEnd, &start, &end, sharedOccCount + ALPHABET_SIZE)) {

				long numAlignNew = *numAlign;
				int tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						bwt, Queries, queryId, 0, sizeX, 1, &start,&end, &numAlignNew);

				if(numAlignNew == -1){
					*numAlign = -1;
					return 0;
				}

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
					*numAlign = numAlignNew;
				}
			}
		}
	}

	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end-start<revEnd-revStart) {
		revStart=revEnd-(end-start);
	} else if (end-start>revEnd-revStart) {
		start=end-(revEnd-revStart);
	}

	if(sizeY - sizeX > TABLE_SIZE) {
		ForwardSearch(revbwt, Queries, sizeX + TABLE_SIZE ,sizeY,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE);
		i = sizeY;

	}
	else
		i = sizeX + TABLE_SIZE;
	for (; i<sizeZ && start<=end; i++)
	{
		c = getCharFromQueries(Queries, i);
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		for (misBase = 0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start <= temp_end)
				ForwardSearch(revbwt, Queries, i+1,sizeZ, &temp_rev_start, &temp_rev_end, 
						&temp_start, &temp_end, sharedOccCount + ALPHABET_SIZE);
			if(temp_start <= temp_end){

				long numAlignNew = *numAlign; 
				int tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						bwt, Queries, queryId,0,sizeX,1,&temp_start, &temp_end, 
						&numAlignNew);

				if(numAlignNew == -1){
					*numAlign = -1;
					return 0;
				}

				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
					*numAlign = numAlignNew;
				}

			}

		}
		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}


	return numOutput;
}



__device__ int BWT2MismatchSearch_CaseD2(uint3* result, int maxNumOutput,
		unsigned int * Queries, unsigned int queryId, BWT *bwt, BWT *revbwt,
		unsigned int *sharedOccCount)
{

	int numOutput = 0;

	unsigned int temp_start=0,temp_end=0,temp_rev_start=0,temp_rev_end=0;
#if (OCC_COUNT_PER_THREAD >= 4)
	unsigned int*occCount = sharedOccCount;
#else
	unsigned int occCount[ALPHABET_SIZE];
#endif

	unsigned int occCount_start[ALPHABET_SIZE];
	unsigned int occCount_end[ALPHABET_SIZE];
	//	unsigned int sizeX = floor((double)QUERY_LENGTH * 0.3f);
	//	unsigned int sizeY = sizeX * 2;
	//	unsigned int sizeZ = QUERY_LENGTH;

	unsigned long long packedPattern = 0;
	//unsigned long long l_packedPattern = 0;
	//unsigned long long r_packedPattern = 0;
	unsigned long long rev_packedPattern = 0;
	//unsigned long long rev_l_packedPattern = 0;
	//unsigned long long rev_r_packedPattern = 0;

	unsigned long long ALLONE = (1<<(bwt->lookupTable->tableSize*2))-1;
	unsigned long long mask;

	unsigned int i;
	char c, misBase;
	unsigned int start, end, revStart, revEnd;

	packedPattern = 0;
	//l_packedPattern = 0;
	//r_packedPattern = 0;
	rev_packedPattern = 0;
	//rev_l_packedPattern = 0;
	//rev_r_packedPattern = 0;

	for (i = 0; i < bwt->lookupTable->tableSize ; i++) {
		packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+i);
		packedPattern |= (c & 3);
	}
	for (i = 0; i <revbwt->lookupTable->tableSize ; i++) {
		rev_packedPattern<<=2;
		c = getCharFromQueries(Queries, sizeX+revbwt->lookupTable->tableSize-i-1);
		rev_packedPattern |= (c & 3);
	}

	//For error happen in lookup range....
	for (i = sizeY - sizeX; i< bwt->lookupTable->tableSize ;i++) {
		c = getCharFromQueries(Queries, sizeX + i);
		for (misBase=0;misBase<4;misBase++) {
			if (misBase == c) continue;
			unsigned long long err_packedPattern = packedPattern;
			unsigned long long err_rev_packedPattern = rev_packedPattern;
			unsigned int bitPos = (bwt->lookupTable->tableSize-i-1)*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_packedPattern |= (3 << bitPos);
			err_packedPattern &= mask;

			bitPos = i*2;
			mask = ALLONE - (3 << bitPos);
			mask |= misBase << bitPos;
			err_rev_packedPattern |= (3 << bitPos);
			err_rev_packedPattern &= mask;

			start = err_packedPattern ? bwt->lookupTable->table[err_packedPattern-1]+1 : 1;
			end = bwt->lookupTable->table[err_packedPattern];
			revStart = err_rev_packedPattern ? revbwt->lookupTable->table[err_rev_packedPattern-1]+1 : 1;
			revEnd = revbwt->lookupTable->table[err_rev_packedPattern];

			//Lookup SA correction...
			if (end - start < revEnd - revStart) {
				revStart = revEnd - ( end - start );
			} else if ( end - start >revEnd - revStart) {
				start = end - (revEnd - revStart);
			}
			if( start <= end && ForwardSearch(revbwt, Queries, sizeX+revbwt->lookupTable->tableSize, QUERY_LENGTH, 
						&revStart, &revEnd, &start, &end, sharedOccCount + ALPHABET_SIZE )) {


				int tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						bwt, Queries, queryId, 0, sizeX, 1, &start,&end);
				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
				}

				//BackwardSearch1Error(bwt, Queries, queryId, 0, sizeX, 1, &start,&end,saInterval);
			}
		}
	}

	start = packedPattern ? bwt->lookupTable->table[packedPattern-1]+1 : 1;
	end = bwt->lookupTable->table[packedPattern];
	revStart = rev_packedPattern ? revbwt->lookupTable->table[rev_packedPattern-1]+1 : 1;
	revEnd = revbwt->lookupTable->table[rev_packedPattern];

	//Lookup SA correction...
	if (end-start<revEnd-revStart) {
		revStart=revEnd-(end-start);
	} else if (end-start>revEnd-revStart) {
		start=end-(revEnd-revStart);
	}

	if(sizeY - sizeX > TABLE_SIZE) {
		ForwardSearch(revbwt, Queries, sizeX + TABLE_SIZE ,sizeY,
				&revStart,&revEnd,&start,&end, sharedOccCount + ALPHABET_SIZE);
		i = sizeY;

	}
	else
		i = sizeX + TABLE_SIZE;
	for (; i<sizeZ && start<=end; i++)
	{
		c = getCharFromQueries(Queries, i);
		bwtCalcAll(revbwt, &revStart, &revEnd,occCount_start, occCount_end, occCount);
		/*
#if LAYER_INDEX_1
GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revStart, revbwt->inverseSa0, occCount_start) ;
GPUBWTAllOccValue(revbwt->bwtCode, revbwt->GPUOccValue, revEnd + 1, revbwt->inverseSa0, occCount_end) ;
#endif

#if LAYER_INDEX_2
BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, revStart,
revbwt->inverseSa0, occCount_start) ;
BWTAllOccValue_2Layer(revbwt->bwtCode, revbwt->occMajorValue, revbwt->occValue, revEnd + 1,
revbwt->inverseSa0, occCount_end) ;
#endif
		//		getAllOccValue(revbwt,revStart, occCount_start);
		//		getAllOccValue(revbwt,revEnd+1, occCount_end);
		occCount[ALPHABET_SIZE-1] = 0;
		for(k=ALPHABET_SIZE-2; k>=0; k--)
		{
		occCount[k] = occCount[k+1] + occCount_end[k+1] - occCount_start[k+1];
		}
		 */
		for (misBase = 0; misBase<ALPHABET_SIZE; misBase++)
		{
			if (misBase == c)
				continue;
			temp_start = start;
			temp_end = end;
			temp_rev_start = revStart;
			temp_rev_end = revEnd;

			bwtCalcWithOccValue(revbwt->countValue,misBase,occCount_start[misBase], occCount_end[misBase], &temp_rev_start, &temp_rev_end);

			temp_end = temp_end - occCount[misBase];
			temp_start = temp_end - (temp_rev_end - temp_rev_start);

			if (temp_start <= temp_end)
				ForwardSearch(revbwt, Queries, i+1,sizeZ, &temp_rev_start, &temp_rev_end, 
						&temp_start, &temp_end, sharedOccCount + ALPHABET_SIZE);
			if(temp_start <= temp_end){

				int tempN = BackwardSearch1Error2(result + numOutput, maxNumOutput - numOutput,
						bwt, Queries, queryId,0,sizeX,1,&temp_start, &temp_end);
				if(tempN == -1) {
					return -1;
				} else {
					numOutput += tempN;
				}

				//BackwardSearch1Error(bwt, Queries, queryId,0,sizeX,1,&temp_start, &temp_end, saInterval);
			}

		}
		bwtCalcWithOccValue(revbwt->countValue,c,occCount_start[c], occCount_end[c], &revStart, &revEnd);

		end = end - occCount[c];
		start = end - (revEnd - revStart);
	}


	return numOutput;
}

/*
__device__ void loadQueries(unsigned int *sharedQueries, unsigned int *Queries)
{
	int iteration,i;
	unsigned int queryId, sharedQueryId;
	for(iteration =0; iteration < ceil(QUERY_NUM / (double)(gridDim.x * blockDim.x)); iteration++) {
		queryId = iteration * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
		if(queryId >= QUERY_NUM)
			break;
		sharedQueryId = iteration * blockDim.x + threadIdx.x;
		for(i=0; i<QUERY_SIZE_IN_WORD; i++)
			sharedQueries[sharedQueryId* QUERY_SIZE_IN_WORD + i] = Queries[queryId * QUERY_SIZE_IN_WORD + i];	
	}
}

__device__ void loadSaInterval(unsigned int *saInterval, unsigned int * d_saInterval)
{
	int iteration,i;
	unsigned int queryId, sharedQueryId;
	for(iteration =0; iteration < ceil(QUERY_NUM / (double)(gridDim.x * blockDim.x)); iteration++) {
		queryId = iteration * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
		if(queryId >= QUERY_NUM)
			break;
		sharedQueryId = iteration * blockDim.x + threadIdx.x;
		for(i=0; i<SA_INTERVAL_PER_QUERY; i++)
			d_saInterval[queryId * SA_INTERVAL_PER_QUERY + i] = saInterval[sharedQueryId* SA_INTERVAL_PER_QUERY + i];	
	}
}

*/
__global__ void CreateOccValue(BWT * bwt, unsigned int *occValueArray)
{
	int i,j,k;
	for(k=0; k< ceil(bwt->textLength / (double)(gridDim.x * blockDim.x * GPU_OCC_INTERVAL)); k++) {
		i = k * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
		if((i*GPU_OCC_INTERVAL) >= bwt->textLength)
			break;
		for(j=0; j<ALPHABET_SIZE; j++) {
#if LAYER_INDEX_1
			occValueArray[i*ALPHABET_SIZE + j] = GPUBWTOccValue(bwt->bwtCode,bwt->GPUOccValue, 
					i*GPU_OCC_INTERVAL, j, bwt->inverseSa0);
#endif

#if LAYER_INDEX_2
			occValueArray[i*ALPHABET_SIZE + j] = BWTOccValue_2Layer(bwt->bwtCode,bwt->occMajorValue,bwt->occValue, 
					i*GPU_OCC_INTERVAL, j, bwt->inverseSa0);


#endif
		}
	}
}


__global__ void exactMatch_kernel(uint3* d_result, unsigned int *d_Queries, const int * d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt )
{
	int numTotalThread = NUM_TOTAL_THREAD;
	unsigned int start, end;
	int i;
	char c;
	unsigned long long packedPattern = 0;
	int queryLen = QUERY_LENGTH;
	unsigned int * query = NULL;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		packedPattern = 0;
		query = d_Queries + QUERY_SIZE_IN_WORD * idx;
		for (i = 0; i < TABLE_SIZE ; i++) {
			packedPattern<<=2;
			c = getCharFromQueries(query, queryLen - TABLE_SIZE + i);
			packedPattern |= (c & 3);
		}
		start = packedPattern ? d_bwt->lookupTable->table[packedPattern-1]+1 : 1;
		end = d_bwt->lookupTable->table[packedPattern];

		if(BackwardSearch(d_bwt, query, 0, QUERY_LENGTH-TABLE_SIZE, &start, &end)){
			d_result[idx].x = idx;
			d_result[idx].y = start;
			d_result[idx].z = end;		
		}
		else 
			d_result[idx].x = 0xFFFFFFFF;
	}
}

__global__ void exactMatch_limit_kernel(long* d_numAlign, uint3* d_result, int* d_numOutput,
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	unsigned int start, end;
	int i;
	char c;
	unsigned long long packedPattern = 0;
	int queryLen = QUERY_LENGTH;
	unsigned int * query = NULL;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		if(d_numAlign[idx] == -1)
			continue;
		packedPattern = 0;
		query = d_Queries + QUERY_SIZE_IN_WORD * idx;
		for (i = 0; i < TABLE_SIZE ; i++) {
			packedPattern<<=2;
			c = getCharFromQueries(query, queryLen - TABLE_SIZE + i);
			packedPattern |= (c & 3);
		}
		start = packedPattern ? d_bwt->lookupTable->table[packedPattern-1]+1 : 1;
		end = d_bwt->lookupTable->table[packedPattern];

		if(BackwardSearch(d_bwt, query, 0, QUERY_LENGTH-TABLE_SIZE, &start, &end)){
			d_result[idx * maxNumOutputPerQuery].x = idx;
			d_result[idx * maxNumOutputPerQuery].y = start;
			d_result[idx * maxNumOutputPerQuery].z = end;
			d_numAlign[idx] += end - start + 1;
			if(d_numAlign[idx] > MAX_NUM_ALIGN_PER_QUERY) {
				d_numAlign[idx] = -1;
				d_numOutput[idx] = 0;	
				continue;
			}
			else {
				d_numOutput[idx] = 1;
			}
		}
		else {
			d_numOutput[idx] = 0;
		}
	}
}

__global__ void BWT2MismatchSearchCaseA_kernel(	uint3* d_result, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryA[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;	
		numOutputCurCase = BWT2MismatchSearch_CaseA2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], 
				&bwt, &revbwt);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		

	}
}


__global__ void BWT2MismatchSearchCaseALimit_kernel(long* d_numAlign, uint3* d_result, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	int queryId = 0;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		queryId = d_QueriesId[idx];
		numOutputCurCase = 0;
		if(d_numAlign[queryId] == -1 || d_numOutput[idx] == -1) {
			continue;
		}

		long numAlignNew = d_numAlign[queryId];
		numOutputCurCase = BWT2MismatchSearch_CaseA2_limit(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, queryId, 
				&bwt, &revbwt, &(numAlignNew));

		if(numAlignNew == -1) {
			d_numAlign[queryId] = -1;
			continue;
		}

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			d_numAlign[queryId] = numAlignNew;
		}

		d_numOutput[idx] += numOutputCurCase;		
	}
}





__global__ void BWT2MismatchSearchCaseB_kernel(	uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];
	//extern __shared__ unsigned int s_queryB[];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryB[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;	
		numOutputCurCase = BWT2MismatchSearch_CaseB2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}




__global__ void BWT2MismatchSearchCaseBLimit_kernel(long* d_numAlign, uint3* d_result, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	int queryId = 0;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		queryId = d_QueriesId[idx];

		if(d_numAlign[queryId] == -1 || d_numOutput[idx] == -1) {
			continue;
		}


		numOutputCurQuery = 0;	
		long numAlignNew = d_numAlign[queryId];
		numOutputCurCase = BWT2MismatchSearch_CaseB2_limit(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, queryId, 
				&bwt, &revbwt,sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD, &(numAlignNew));

		if(numAlignNew == -1) {
			d_numAlign[queryId] = -1;
			continue;
		}

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
			d_numAlign[queryId] = numAlignNew;
		}

		d_numOutput[idx] += numOutputCurCase;		
	}
}




__global__ void BWT2MismatchSearchCaseC_kernel(	uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryC[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;

		numOutputCurCase = BWT2MismatchSearch_CaseC2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt, 
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}


__global__ void BWT2MismatchSearchCaseCLimit_kernel(long* d_numAlign, uint3* d_result, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurCase;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	int queryId = 0;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		queryId = d_QueriesId[idx];

		if(d_numAlign[idx] == -1 || d_numOutput[idx] == -1) {
			continue;
		}


		long numAlignNew = d_numAlign[idx];
		numOutputCurCase = BWT2MismatchSearch_CaseC2_limit(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, queryId, 
				&bwt, &revbwt, sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD,&(numAlignNew));

		if(numAlignNew == -1) {
			d_numAlign[idx] = -1;
			continue;
		}

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			d_numAlign[idx] = numAlignNew;
		}

		d_numOutput[idx] += numOutputCurCase;		
	}
}






__global__ void BWT2MismatchSearchCaseD_kernel(	uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryD[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;
		numOutputCurCase = BWT2MismatchSearch_CaseD2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;		
	}
}

__global__ void BWT2MismatchSearch_4Case_kernel(uint3* d_result, int* d_numOutput, int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	__shared__ unsigned int sharedOccCount[OCC_COUNT_PER_THREAD * MAX_THREAD_NUM];

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		//load the query into the shared memory
		/*for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
		  s_queryD[threadIdx.x*QUERY_SIZE_IN_WORD + wordId] = 
		  d_Queries[idx*QUERY_SIZE_IN_WORD + wordId];
		  }*/

		numOutputCurQuery = 0;
		numOutputCurCase = 0;
		numOutputCurCase = BWT2MismatchSearch_CaseA2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseB2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseC2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);
		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}

		numOutputCurCase = BWT2MismatchSearch_CaseD2(d_result + maxNumOutputPerQuery*idx + numOutputCurQuery,
				maxNumOutputPerQuery - numOutputCurQuery,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx] = numOutputCurQuery;
	}
}

__global__ void BWT2MismatchSearchCaseDLimit_kernel(long* d_numAlign, uint3* d_result, int* d_numOutput, 
		int maxNumOutputPerQuery,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	int queryId = 0;
	__shared__ unsigned int sharedOccCount[MAX_THREAD_NUM * OCC_COUNT_PER_THREAD];
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		queryId = d_QueriesId[idx];

		if(d_numAlign[queryId] == -1 || d_numOutput[idx] == -1) {
			continue;
		}


		numOutputCurQuery = 0;	
		long numAlignNew = d_numAlign[queryId];
		numOutputCurCase = BWT2MismatchSearch_CaseD2_limit(d_result + maxNumOutputPerQuery*idx + d_numOutput[idx],
				maxNumOutputPerQuery - d_numOutput[idx],
				d_Queries + QUERY_SIZE_IN_WORD*idx, queryId, 
				&bwt, &revbwt, sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD, &(numAlignNew));

		if(numAlignNew == -1) {
			d_numAlign[queryId] = -1;
			continue;
		}

		if(numOutputCurCase == -1) {
			d_numOutput[idx] = -1;
			continue; 
		} else {
			numOutputCurQuery += numOutputCurCase;
			d_numAlign[queryId] = numAlignNew;
		}

		d_numOutput[idx] += numOutputCurCase;		
	}
}




__global__ void BWT1MismatchForward_shareBuf_kernel(uint3* d_result, int2* d_numOutput, int maxNumOutputPerThread,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{
	extern __shared__ unsigned int sharedOccCount[];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;
		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurCase = BWT1MismatchSearch_Forward(d_result + maxNumOutputPerThread* threadId + resultCount,
				maxNumOutputPerThread - resultCount,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt,
				sharedOccCount + OCC_COUNT_PER_THREAD * threadIdx.x);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue; 
		} else{
			numOutputCurQuery += numOutputCurCase;
			d_numOutput[idx].x = resultCount;          
			d_numOutput[idx].y = resultCount + numOutputCurCase;
			resultCount += numOutputCurCase;
		}

	}
}


__global__ void BWT1MismatchBackward_shareBuf_kernel(uint3* d_result, int2* d_numOutput, int maxNumOutputPerThread,
		unsigned int * Queries, unsigned int queryNum, BWT *bwt, BWT *revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT * ptr_bwt = bwt;
	BWT * ptr_revbwt = revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {
		numOutputCurQuery = 0;
		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurCase = BWT1MismatchSearch_Backward(d_result + maxNumOutputPerThread* threadId + resultCount,
				maxNumOutputPerThread - resultCount,
				Queries + QUERY_SIZE_IN_WORD*idx, idx, ptr_bwt, ptr_revbwt);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD  > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue; 
		} else{
			numOutputCurQuery += numOutputCurCase;
			d_numOutput[idx].x = resultCount;          
			d_numOutput[idx].y = resultCount + numOutputCurCase;
			resultCount += numOutputCurCase;
		}

	}
}


__global__ void BWT2MismatchSearchCaseA_shareBuf_kernel(        uint3* d_result, int2* d_numOutput,
		int maxNumOutputPerThread,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;
	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurQuery = 0;  

		numOutputCurCase = BWT2MismatchSearch_CaseA2(d_result + threadId * maxNumOutputPerThread + resultCount,
				maxNumOutputPerThread - resultCount,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], 
				&bwt, &revbwt);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD  > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue; 
		} else{
			numOutputCurQuery += numOutputCurCase;
			d_numOutput[idx].x = resultCount;          
			d_numOutput[idx].y = resultCount + numOutputCurCase;
			resultCount += numOutputCurCase;
		}
	}
}

__global__ void BWT2MismatchSearchCaseB_shareBuf_kernel(        uint3* d_result, int2* d_numOutput, int maxNumOutputPerThread,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	extern __shared__ unsigned int sharedOccCount[];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		numOutputCurQuery = 0;
		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurCase = BWT2MismatchSearch_CaseB2(d_result + maxNumOutputPerThread * threadId + resultCount,
				maxNumOutputPerThread - resultCount,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt,
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD  > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;

			d_numOutput[idx].x = resultCount;          
			d_numOutput[idx].y = resultCount + numOutputCurCase;
			resultCount += numOutputCurCase;
		}
	}
}

__global__ void BWT2MismatchSearchCaseC_shareBuf_kernel(        uint3* d_result, int2* d_numOutput, int maxNumOutputPerThread,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	extern __shared__ unsigned int sharedOccCount[];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		numOutputCurQuery = 0;
		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurCase = BWT2MismatchSearch_CaseC2(d_result + maxNumOutputPerThread * threadId + resultCount,
				maxNumOutputPerThread - resultCount,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt, 
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD  > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx].x = resultCount;          
		d_numOutput[idx].y = resultCount + numOutputCurCase;
		resultCount += numOutputCurCase;
	}
}

__global__ void BWT2MismatchSearchCaseD_shareBuf_kernel(        uint3* d_result, int2* d_numOutput, int maxNumOutputPerThread,
		unsigned int * d_Queries, const int* d_QueriesId,
		unsigned int queryNum, BWT *d_bwt, BWT *d_revbwt)
{
	extern __shared__ unsigned int sharedOccCount[];
	int numTotalThread = NUM_TOTAL_THREAD;
	int numOutputCurQuery, numOutputCurCase;
	BWT bwt = *d_bwt;
	BWT revbwt = *d_revbwt;
	bool bufferFull = false;
	int resultCount = 0;
	int threadId = blockIdx.x * blockDim.x + threadIdx.x;

	for(int idx = GLOBAL_THREAD_OFFSET; idx < queryNum; idx += numTotalThread) {

		numOutputCurQuery = 0;
		if(bufferFull){
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		}
		numOutputCurCase = BWT2MismatchSearch_CaseD2(d_result + maxNumOutputPerThread * threadId + resultCount,
				maxNumOutputPerThread - resultCount,
				d_Queries + QUERY_SIZE_IN_WORD*idx, d_QueriesId[idx], &bwt, &revbwt, 
				sharedOccCount + threadIdx.x * OCC_COUNT_PER_THREAD);

		if(numOutputCurCase == -1) {
			if(resultCount + RESULT_BUFFER_THRESHOLD  > maxNumOutputPerThread)
				bufferFull = true;
			d_numOutput[idx].x = -1;
			d_numOutput[idx].y = -1;
			continue;
		} else {
			numOutputCurQuery += numOutputCurCase;
		}
		d_numOutput[idx].x = resultCount;          
		d_numOutput[idx].y = resultCount + numOutputCurCase;
		resultCount += numOutputCurCase;
	}

}

int getNumOutputPerThread(const uint64 resultBufSize, const int numBlock, const int numThread, const int numQuery)
{
	int threadNum = numThread * numBlock;
	uint64 bufSizePerThread = floor(resultBufSize / (double)threadNum);
	int numOutputPerThread = floor(bufSizePerThread / (double)sizeof(uint3));
	return numOutputPerThread;
}


int getNumOutputPerQuery_shareBuf(const uint64 resultBufSize, const int queryNum, const int numBlock, const int numThread) {

	int queryNum0 = ceil(queryNum / (double) (numBlock * numThread)) * numBlock * numThread;
	int bufSizePerQuery = floor(resultBufSize/(double)queryNum0);
	int numOutputPerQuery = floor(bufSizePerQuery/(double)sizeof(uint3));

	return numOutputPerQuery;
}


int getNumOutputPerQuery(const uint64 resultBufSize, const int queryNum) {

	int bufSizePerQuery = floor(resultBufSize/(double)queryNum);
	int numOutputPerQuery = floor(bufSizePerQuery/(double)sizeof(uint3));

	return numOutputPerQuery;
}


__global__
void setLabel_kernel(int* d_temp, const int* d_numResultPerQuery, const int numQuery) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numResultPerQuery[idx] == -1) {
			d_temp[idx] = 1;
		}
	}
}


__global__
void setLabelLimit_kernel(const long* d_numAlign, const int* d_queryId,
		int* d_temp, const int* d_numResultPerQuery, const int numQuery) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numAlign[d_queryId[idx]] != -1 && d_numResultPerQuery[idx] == -1) {
			d_temp[idx] = 1;
		}
	}
}



__global__ void collectUndoQuery_kernel(unsigned int* d_undoQuery, int* d_undoQueryId, 
		const unsigned int* d_query, const int* d_queryId, 
		const int numQuery,
		const int* d_numResultPerQuery, const int* d_pos) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numResultPerQuery[idx] == -1) {
			const int pos = d_pos[idx];
			for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
				d_undoQuery[pos*QUERY_SIZE_IN_WORD + wordId] = 
					d_query[idx*QUERY_SIZE_IN_WORD + wordId];
			}
			d_undoQueryId[pos] = d_queryId[idx];
		}
	}
}

__global__
void setLabel_shareBuf_kernel(int* d_temp, const int2* d_numResultPerQuery, const int numQuery) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numResultPerQuery[idx].x == -1) {
			d_temp[idx] = 1;
		}
	}
}



__global__ void collectUndoQuery_shareBuf_kernel(unsigned int* d_undoQuery, int* d_undoQueryId,
		const unsigned int* d_query, const int* d_queryId,
		const int numQuery,
		const int2* d_numResultPerQuery, const int* d_pos) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numResultPerQuery[idx].x == -1) {
			const int pos = d_pos[idx];
			for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
				d_undoQuery[pos*QUERY_SIZE_IN_WORD + wordId] =
					d_query[idx*QUERY_SIZE_IN_WORD + wordId];
			}
			d_undoQueryId[pos] = d_queryId[idx];
		}
	}
}


__global__ void collectUndoQueryLimit_kernel(const long* d_numAlign, unsigned int* d_undoQuery, int* d_undoQueryId, 
		const unsigned int* d_query, const int* d_queryId, 
		const int numQuery,
		const int* d_numResultPerQuery, const int* d_pos) {
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		if(d_numAlign[d_queryId[idx]] != -1 && d_numResultPerQuery[idx] == -1) {
			const int pos = d_pos[idx];
			for(int wordId = 0; wordId < QUERY_SIZE_IN_WORD; wordId++) {
				d_undoQuery[pos*QUERY_SIZE_IN_WORD + wordId] = 
					d_query[idx*QUERY_SIZE_IN_WORD + wordId];
			}
			d_undoQueryId[pos] = d_queryId[idx];
		}
	}
}


__global__
void countAlign_shareBuf_kernel(	
		int* d_temp, const uint3* d_result, const int2* d_numResultPerQuery, 
		const int maxNumResultPerThread, const int numQuery, 
		const int* d_queryId, int* d_queryLabel) {

	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;
	int threadId = blockDim.x * blockIdx.x + threadIdx.x;
	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		const int numResult = d_numResultPerQuery[idx].y - d_numResultPerQuery[idx].x;
		if(numResult > 0) {
			int numAlign = 0;
			d_queryLabel[d_queryId[idx]] = 1;

			for(int i = d_numResultPerQuery[idx].x; i < d_numResultPerQuery[idx].y; i++) {
				numAlign += (d_result[threadId * maxNumResultPerThread + i].z 
						- d_result[threadId * maxNumResultPerThread + i].y + 1);
			}
			//		if(numAlign < 0)
			//			printf("idx: %d\t numAlign:%d\n",idx, numAlign);
			d_temp[idx] = numAlign;
		}
	}
}


__global__
void countAlign_kernel(	int* d_temp, const uint3* d_result, const int* d_numResultPerQuery, 
		const int maxNumResultPerQuery, const int numQuery, 
		const int* d_queryId, int* d_queryLabel) {

	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	for(int idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) {
		const int numResult = d_numResultPerQuery[idx];
		if(numResult > 0) {
			int numAlign = 0;
			d_queryLabel[d_queryId[idx]] = 1;

			for(int i = 0; i < numResult; i++) {
				numAlign += (d_result[idx*maxNumResultPerQuery + i].z - 
						d_result[idx*maxNumResultPerQuery + i].y + 1);
			}
			d_temp[idx] = numAlign;
		}
	}
}

void countAlignGPU_shareBuf(	uint64* numAlignment, 
		const uint3* d_result, const int2* d_numResultPerQuery, 
		const int maxNumResultPerThread, const int numQuery,
		int* d_temp, 
		const int* d_queryId, int* d_queryLabel,int numBlock, int numThread) {
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	printf("count align gpu, maxNumResultPerThread: %d, numQuery: %d\n", maxNumResultPerThread, numQuery);
	countAlign_shareBuf_kernel<<<numBlock, numThread>>>(d_temp, d_result, d_numResultPerQuery, 
			maxNumResultPerThread, numQuery, 
			d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
}



//do not collect for the next pass, just count the alignment
void countAlignGPU(	uint64* numAlignment, 
		const uint3* d_result, const int* d_numResultPerQuery, 
		const int maxNumResultPerQuery, const int numQuery,
		int* d_temp, 
		const int* d_queryId, int* d_queryLabel) {
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	countAlign_kernel<<<1024, 512>>>(d_temp, d_result, d_numResultPerQuery, 
			maxNumResultPerQuery, numQuery, 
			d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
}


int getNumMatchedQuery(int* d_queryLabel, const int numQuery) {
	thrust::device_ptr<int> dptr_queryLabel(d_queryLabel);
	return thrust::reduce(dptr_queryLabel, dptr_queryLabel + numQuery);
}


void collectUndoneGPU_shareBuf(	uint64* numAlignment, 
		unsigned int* d_undoQuery, int* d_undoQueryId, int* numUndo, 
		const uint3* d_result, const int2* d_numResultPerQuery, const int maxNumResultPerThread, 
		const unsigned int* d_query, const int* d_queryId, const int numQuery,
		int* d_temp /*#=numQuery*/, int* d_queryLabel, int numBlock, int numThread) {
	(*numAlignment) = 0;
	(*numUndo) = 0;

	//count the #query should be collected for the next pass
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));
	setLabel_shareBuf_kernel<<<numBlock, numThread>>>(d_temp, d_numResultPerQuery, numQuery);
	cutilCheckMsg("setLabel_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);

	int* h_temp = (int*)malloc(sizeof(int)*numQuery);
	FROMGPU(h_temp, d_temp, sizeof(int)*numQuery);
#if __DEBUG__
	for(int i = 0; i < numQuery; i++) {
		if(h_temp[i] != 0 && h_temp[i] != 1) {
			printf("%d, %d\n", i, h_temp[i]);
		}
	}
#endif
	(*numUndo) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	thrust::exclusive_scan(dptr_temp, dptr_temp + numQuery, dptr_temp);

	//collect the query for the next pass
	collectUndoQuery_shareBuf_kernel<<<numBlock, numThread>>>(d_undoQuery, d_undoQueryId, 
			d_query, d_queryId, numQuery,
			d_numResultPerQuery, d_temp);
	cutilCheckMsg("collectUndoQuery_kernel");

	//count the alignment
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	countAlign_shareBuf_kernel<<<numBlock, numThread>>>(d_temp, d_result, d_numResultPerQuery, 
			maxNumResultPerThread, numQuery, 
			d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);

}


void collectUndoneGPU(	uint64* numAlignment, 
		unsigned int* d_undoQuery, int* d_undoQueryId, int* numUndo, 
		const uint3* d_result, const int* d_numResultPerQuery, const int maxNumResultPerQuery, 
		const unsigned int* d_query, const int* d_queryId, const int numQuery,
		int* d_temp /*#=numQuery*/, int* d_queryLabel) {
	(*numAlignment) = 0;
	(*numUndo) = 0;

	//count the #query should be collected for the next pass
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));
	setLabel_kernel<<<1024, 512>>>(d_temp, d_numResultPerQuery, numQuery);
	cutilCheckMsg("setLabel_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	(*numUndo) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	thrust::exclusive_scan(dptr_temp, dptr_temp + numQuery, dptr_temp);

	//collect the query for the next pass
	collectUndoQuery_kernel<<<1024, 512>>>(d_undoQuery, d_undoQueryId, 
			d_query, d_queryId, numQuery,
			d_numResultPerQuery, d_temp);
	cutilCheckMsg("collectUndoQuery_kernel");

	//count the alignment
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	countAlign_kernel<<<1024, 512>>>(d_temp, d_result, d_numResultPerQuery, 
			maxNumResultPerQuery, numQuery, 
			d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
}





void collectUndoneGPULimit(long* d_numAlign, 
		unsigned int* d_undoQuery, int* d_undoQueryId, int* numUndo, 
		const uint3* d_result, const int* d_numResultPerQuery, const int maxNumResultPerQuery, 
		const unsigned int* d_query, const int* d_queryId, const int numQuery,
		int* d_queryLabel, int* d_temp) {
	(*numUndo) = 0;

	//count the #query should be collected for the next pass
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));
	setLabelLimit_kernel<<<1024, 512>>>(d_numAlign, d_queryId, d_temp, d_numResultPerQuery, numQuery);
	cutilCheckMsg("setLabel_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	(*numUndo) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	thrust::exclusive_scan(dptr_temp, dptr_temp + numQuery, dptr_temp);

	//collect the query for the next pass
	collectUndoQueryLimit_kernel<<<1024, 512>>>(d_numAlign, d_undoQuery, d_undoQueryId, 
			d_query, d_queryId, numQuery,
			d_numResultPerQuery, d_temp);
	cutilCheckMsg("collectUndoQuery_kernel");
	/*	
	//count the alignment
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	countAlign_kernel<<<1024, 512>>>(d_temp, d_result, d_numResultPerQuery, 
	maxNumResultPerQuery, numQuery, 
	d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	 */
}



void collectUndone(	uint64* numAlignment, 
		unsigned int* undoQuery, int* undoQueryId,
		int* numUndo,
		const uint3* result, const int* numResultPerQuery, const int maxNumResultPerQuery,
		const unsigned int* query, const int* inQueryId,
		const int numQuery) {
	(*numAlignment) = 0;
	(*numUndo) = 0;

	for(int idx = 0; idx < numQuery; idx++) {

		int queryId = inQueryId[idx];

		//should be processed in the next pass
		if(numResultPerQuery[idx] == -1) {
			memcpy(undoQuery + QUERY_SIZE_IN_WORD*(*numUndo), query + QUERY_SIZE_IN_WORD*idx, 
					sizeof(unsigned int)*QUERY_SIZE_IN_WORD);
			undoQueryId[*numUndo]  = queryId;

			(*numUndo)++;
		} else {
			//calculate the number of alignment for the done queries
			assert(numResultPerQuery[idx] >= 0 && numResultPerQuery[idx] <= maxNumResultPerQuery);

			for(int i = 0; i < numResultPerQuery[idx]; i++) {
				(*numAlignment) += (result[idx*maxNumResultPerQuery + i].z - 
						result[idx*maxNumResultPerQuery + i].y + 1);
			}
		}
	}
}
unsigned int * SortQueriesBackward(const unsigned int * Queries, const unsigned int queryNum)
{
	int i,j;
	size_t N = 1000000;
	cudaEvent_t start;
	cudaEvent_t end;
	float elapsed_time = 0;
	unsigned int * sortedQuery = (unsigned int *) malloc(sizeof(int) * QUERY_SIZE_IN_WORD * queryNum);

	cudaEventCreate(&start);
	cudaEventCreate(&end);

	thrust::host_vector<unsigned int> host_queryId(queryNum);
	thrust::host_vector<unsigned int> host_key(queryNum);
	for(i=0; i<queryNum; i++){
		host_queryId[i] = i;
		host_key[i] = Queries[i*QUERY_SIZE_IN_WORD + QUERY_SIZE_IN_WORD - 1];
	}
	thrust::device_vector<unsigned int> device_queryId = host_queryId;
	thrust::device_vector<unsigned int> device_key = host_key;


	cudaEventRecord(start,0);
	thrust::sort_by_key(device_key.begin(), device_key.end(), device_queryId.begin());
	cudaThreadSynchronize();

	cudaEventRecord(end,0);
	cudaEventSynchronize(end);

	thrust::copy(device_queryId.begin(), device_queryId.end(), host_queryId.begin());

	std::cout << "Queries sort took " << elapsed_time * N << " seconds" << std::endl;
	for(i=0; i < queryNum; i++ ){
		for(j=0; j < QUERY_SIZE_IN_WORD; j++) {
			sortedQuery[i * QUERY_SIZE_IN_WORD + j] = Queries[host_queryId[i] * QUERY_SIZE_IN_WORD + j];
		}
	}

	return sortedQuery;
}
unsigned int * SortQueriesForward(const unsigned int * Queries,const unsigned int queryNum)
{
	int i,j;
	size_t N = 1000000;
	cudaEvent_t start;
	cudaEvent_t end;
	float elapsed_time = 0;
	unsigned int * sortedQuery = (unsigned int *) malloc(sizeof(int) * QUERY_SIZE_IN_WORD * queryNum);

	cudaEventCreate(&start);
	cudaEventCreate(&end);

	thrust::host_vector<unsigned int> host_queryId(queryNum);
	thrust::host_vector<unsigned int> host_key(queryNum);
	for(i=0; i<queryNum; i++){
		host_queryId[i] = i;
		host_key[i] = Queries[i*QUERY_SIZE_IN_WORD];
	}
	thrust::device_vector<unsigned int> device_queryId = host_queryId;
	thrust::device_vector<unsigned int> device_key = host_key;


	cudaEventRecord(start,0);
	thrust::sort_by_key(device_key.begin(), device_key.end(), device_queryId.begin());
	cudaThreadSynchronize();

	cudaEventRecord(end,0);
	cudaEventSynchronize(end);

	thrust::copy(device_queryId.begin(), device_queryId.end(), host_queryId.begin());

	std::cout << "Queries sort took " << elapsed_time * N << " seconds" << std::endl;

	for(i=0; i < queryNum; i++ ){
		for(j=0; j < QUERY_SIZE_IN_WORD; j++) {
			sortedQuery[i * QUERY_SIZE_IN_WORD + j] = Queries[host_queryId[i] * QUERY_SIZE_IN_WORD + j];
		}
	}

	return sortedQuery;
}

__global__ void getResultNumPerThread_kernel(const int2* d_numResultPerQuery, uint32 numQuery, int* d_numResultPerThread)
{
	const int numTotalThread = NUM_TOTAL_THREAD;
	const int globalThreadOffset = GLOBAL_THREAD_OFFSET;

	int idx;
	d_numResultPerThread[globalThreadOffset] = 0;
	for(idx = globalThreadOffset; idx < numQuery; idx += numTotalThread) 
	{
		if(d_numResultPerThread[globalThreadOffset] == d_numResultPerQuery[idx].x)
			d_numResultPerThread[globalThreadOffset] = d_numResultPerQuery[idx].y;
	}
	if(d_numResultPerThread[globalThreadOffset] == -1)
		d_numResultPerThread[globalThreadOffset] = 0;
}

void collectUndoneGPU_shareBuf(	uint64* numAlignment, 
		unsigned int* d_undoQuery, int* d_undoQueryId, int* numUndo, 
		const uint3* d_result, const int2* d_numResultPerQuery, const int maxNumResultPerThread, 
		const unsigned int* d_query, const int* d_queryId, const int numQuery,
		int* d_temp /*#=numQuery*/, int* d_queryLabel, int numBlock, int numThread, 
		uint3* h_result, uint32* h_resultNum) {
	(*numAlignment) = 0;
	(*numUndo) = 0;

	//copy and compact the result on the CPU. TODO: compaction on the GPU
	cutilSafeCall(cudaThreadSynchronize());
	unsigned int timer = 0;
	startTimer(&timer);
	uint3* h_resultStart = h_result + (*h_resultNum);
	int* d_numResultPerThread;
	uint3* d_collectedResult;
	cutilSafeCall(cudaMalloc(&d_numResultPerThread, sizeof(int)*numBlock*numThread));
	cutilSafeCall(cudaMalloc(&d_collectedResult, sizeof(uint3)*30000000));
	cutilSafeCall(cudaMemset(d_collectedResult, 0, sizeof(uint3)*30000000));
	int numTotalResult = 0;
	printf("buffer size: %u\n",  sizeof(int3) *numBlock*numThread);

	
	getResultNumPerThread_kernel<<<numBlock, numThread>>>(d_numResultPerQuery, numQuery, d_numResultPerThread);

	unsigned int timer1 = 0;
	startTimer(&timer1);
	numTotalResult = compaction(d_collectedResult, d_result, numBlock*numThread, maxNumResultPerThread, d_numResultPerThread);
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer1, "compaction kernel");	

	FROMGPU(h_resultStart, d_collectedResult, sizeof(uint3)*numTotalResult);

	cutilSafeCall(cudaFree(d_numResultPerThread));
	cutilSafeCall(cudaFree(d_collectedResult));
	(*h_resultNum) += numTotalResult;
	cutilSafeCall(cudaThreadSynchronize());
	printf("numTotalResult = %d\n", numTotalResult);
//	printf("numTotalResult_cpu = %d\n", numTotalResult_cpu);
	endTimer(&timer, "result compaction");

	//count the #query should be collected for the next pass
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));
	setLabel_shareBuf_kernel<<<numBlock, numThread>>>(d_temp, d_numResultPerQuery, numQuery);
	cutilCheckMsg("setLabel_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	int* h_temp = (int*)malloc(sizeof(int)*numQuery);
	FROMGPU(h_temp, d_temp, sizeof(int)*numQuery);
	(*numUndo) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	thrust::exclusive_scan(dptr_temp, dptr_temp + numQuery, dptr_temp);

	//collect the query for the next pass
	collectUndoQuery_shareBuf_kernel<<<numBlock, numThread>>>(d_undoQuery, d_undoQueryId, 
			d_query, d_queryId, numQuery,
			d_numResultPerQuery, d_temp);
	cutilCheckMsg("collectUndoQuery_kernel");

	//count the alignment
	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));	
	countAlign_shareBuf_kernel<<<numBlock, numThread>>>(d_temp, d_result, d_numResultPerQuery, 
			maxNumResultPerThread, numQuery, 
			d_queryId, d_queryLabel);
	cutilCheckMsg("countAlign_kernel");
	(*numAlignment) = thrust::reduce(dptr_temp, dptr_temp + numQuery);

}

__global__ void getNumResultPerQuery_kernel(int* numResultPerQuery,long* numAlign,int maxNumResultPerQuery, uint32 numQuery) 
{
	uint32 idx = blockDim.x * blockIdx.x + threadIdx.x;
	while(idx < numQuery) {
		if(numResultPerQuery[idx] == -1 || numAlign[idx] == -1 || numResultPerQuery[idx] >= maxNumResultPerQuery)
			numResultPerQuery[idx] = 0;
		idx += blockDim.x * gridDim.x ;
	}
}

void collectUndoneGPULimit(long* d_numAlign,  
		const uint3* d_result, int* d_numResultPerQuery, const int maxNumResultPerQuery, 
		const unsigned int* d_query, const int* d_queryId, const int numQuery,
		uint3* h_result, uint32* h_resultNum) {

	uint3* h_resultStart = h_result + (*h_resultNum);
	uint3* d_collectedResult;
	cutilSafeCall(cudaMalloc(&d_collectedResult, sizeof(uint3)*30000000));
	cutilSafeCall(cudaMemset(d_collectedResult, 0, sizeof(uint3)*30000000));
	int numTotalResult = 0;
	int numElementCopy = 0;

	uint32 timer = 0;
/*	
	long* h_numAlign = (long*)malloc(numQuery * sizeof(long));
	uint32* h_numPerQuery = (uint32*)malloc(numQuery * sizeof(int));
	FROMGPU(h_numAlign, d_numAlign, sizeof(long)*numQuery);
	FROMGPU(h_numPerQuery, d_numResultPerQuery, sizeof(int) * numQuery);

	int tempCount = 0;
	for(int i=0; i < numQuery; i++) {
		if(h_numAlign[i] >0 && h_numPerQuery[i] > 0)
			tempCount += h_numPerQuery[i];
	}
	printf("sa count: %lu\n", tempCount);

	free(h_numAlign);
	free(h_numPerQuery);
*/
	startTimer(&timer);
	getNumResultPerQuery_kernel<<<128,64>>>(d_numResultPerQuery, d_numAlign, maxNumResultPerQuery, numQuery); 

	printf("maxNumResultPerQuery: %u\n", maxNumResultPerQuery);
	printf("numQuery = %d\n", numQuery);

	unsigned int timer1 = 0;
	startTimer(&timer1);	
	printf("2. numQuery = %d\n", numQuery);
	numTotalResult = compaction(d_collectedResult, d_result, numQuery, maxNumResultPerQuery, d_numResultPerQuery);
	cutilSafeCall(cudaThreadSynchronize());
	endTimer(&timer1, "compaction kernel");
	(*h_resultNum) += numTotalResult;
	
	FROMGPU(h_resultStart, d_collectedResult, sizeof(uint3)*numTotalResult);

	cutilSafeCall(cudaFree(d_collectedResult));
	cutilSafeCall(cudaThreadSynchronize());

	printf("numTotalResult = %d\n", numTotalResult);
	endTimer(&timer, "result compaction");
/*
	//count the #query should be collected for the next pass

	cutilSafeCall(cudaMemset(d_temp, 0, sizeof(int)*numQuery));
	setLabelLimit_kernel<<<1024, 512>>>(d_numAlign, d_queryId, d_temp, d_numResultPerQuery, numQuery);
	cutilCheckMsg("setLabel_kernel");
	thrust::device_ptr<int> dptr_temp(d_temp);
	(*numUndo) = thrust::reduce(dptr_temp, dptr_temp + numQuery);
	thrust::exclusive_scan(dptr_temp, dptr_temp + numQuery, dptr_temp);

	//collect the query for the next pass
	collectUndoQueryLimit_kernel<<<1024, 512>>>(d_numAlign, d_undoQuery, d_undoQueryId, 
			d_query, d_queryId, numQuery,
			d_numResultPerQuery, d_temp);
	cutilCheckMsg("collectUndoQuery_kernel");
*/
}




void collectExactResult(uint3* d_result, uint32 numQuery, uint3* h_result, uint32* h_resultNum) {

	assert(*h_resultNum == 0);
	uint3* h_resultStart = h_result + (*h_resultNum);

	int n = 0;
	FROMGPU(h_resultStart, d_result, sizeof(uint3)*numQuery);
	int numCollect = 0;
	for(uint32 i = 0; i < numQuery; i++) {
		if(h_resultStart[i].x < numQuery) {
			h_resultStart[numCollect] = h_resultStart[i];
			n += h_resultStart[i].z - h_resultStart[i].y + 1;
			assert(numCollect <= i);
			numCollect += 1;
		}
	}
	(*h_resultNum) += numCollect;

	printf("======> n = %d\n", n);

}
